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C/^ . Abstract. A theoretical framework is developed for a precise control of spatial 

I patterns in oscillatory media using nonlinear global feedback, where a proper form 

^ ' of the feedback function corresponding to a specific pattern is predicted through the 

analysis of a phase diffusion equation with global coupling. In particular, feedback 
functions that generate the following spatial patterns are analytically given: i) 2- 
cluster states with an arbitrary population ratio, ii) equally populated multi-cluster 
states, and iii) a desynchronized state. Our method is demonstrated numerically by 
' using the Brusselator model in the oscillatory regime. Experimental realization is also 

CN ■ discussed. 
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1. Introduction 



Feedback control is a powerful method of regulating spatio-temporal dynamics and has 
5^ ! been studied in a wide variety of fields including physics, chemistry, biology, and medical 
science [1]. For example, formation of various clustering patterns has been realized in 
the Belousov-Zhabotinsky reaction [2], [3] and in the catalytic CO oxidation reaction 
on Pt [H El E]. The catalytic CO oxidation systems have also been studied for the 
suppression of chemical turbulence [5l [7]. Moreover, considerable attention has been 
paid to feedback devices that suppress the pathological synchronization in the brain of 
Parkinson's disease patients [H El [lOl [IH IHl [13] . 

In many cases systems to be controlled are spatially extended, and reaction-diffusion 
systems provide a good model for the study of pattern controlling. Theoretical analyses 
based on reaction-diffusion systems have been done for the Belousov-Zhabotinsky 
reaction [31 [HI [15] and CO oxidation [161 [HI [18] . However, so far, only empirical control 
has been achieved for such spatially-extended systems, including above-mentioned 
pioneering experimental works [21 [3l [6]; there has been no general theory that 
quantitatively relates feedback inputs to spatial patterns. 
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On the other hand, for discrete oscillator systems, such a quantitative feedback 
control methodology has been established very recently by Kiss, Kori, Hudson, and 
Rusin [ini [20]. Their method is based on a phase model described by 



where (pi (0 < 0j < 2n) is the phase of the oscillator i {i = 1, . . . , N), u is the natural 
frequency, K is the coupling strength, and is called the coupling function. Their 
method utilizes the following facts: the coupling function determines the entire collective 
behavior of the phase model, and any coupling function can be designed by applying 
an appropriately constructed feedback signal to a population of oscillators. Hence 
the population of oscillators can be steered to a desired synchronization behavior by 
taking the following two steps: (i) find a coupling function that results in a desired 
synchronization behavior in (1), and (ii) construct an appropriate feedback signal that 
yields the coupling function. A major advantage of their methodology is that the phase 
model can be constructed from experimentally measurable quantities only; detailed 
information on the intrinsic dynamics of the system is not necessary. Validity and 
robustness of their methodology have been confirmed both experimentally by using 
electro-chemical oscillators [T9l [20] and numerically [20] . 

In this paper, by utilizing the above methodology by Kiss, Kori, Hudson, and 
Rusin, we develop a general theory for the global feedback control of spatially extended 
oscillatory media. Our approach is also based on a phase model. Since the existence 
of diffusive coupling plays a crucial role on the development of spatial patterns in 
oscillatory media, our phase model inevitably includes both diffusive and global coupling, 
in contrast to discrete oscillators. Studying such a phase model, we find coupling 
functions leading to the following spatial patterns characterized by the distribution 
of phases: (i) 2-cluster states with specified population ratios, (ii) equally populated 
multi-cluster states, and (iii) a desynchronized state. Moreover, we propose a new 
nonlinear feedback function without time delay, which is more convenient to design 
various coupling functions than that used in the previous work [19], [20] . We numerically 
demonstrate our proposed method by using a particular reaction-diffusion model and 
reproduce all the above three patterns with theoretically predicted feedback parameters. 

This paper is organized as follows: In section [21 we present the basic idea of our 
control methodology for oscillatory media in detail. In section [3l we give a detailed 
analysis of the phase diffusion equation with special coupling functions that yield the 
above-mentioned three spatial patterns. Numerical demonstration of the theory by using 
the Brusselator is given in section [H Experimental realization is discussed in section [5] 

2. General control methodology 

Our approach to the control of oscillatory media is closely related to the method recently 
proposed for the population of oscillators [I9l[20]. Dynamics of discrete, identical limit- 




(1) 
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cycle oscillators under global feedback is described by the following nonlinear dynamical 
equations: 

^ = F(^..) + |ef^/^K.), (2) 

i=i 

where Ui is the state vector of the i^^ oscillator (z = 1, ■ ■ ■ , A^), is a nonlinear function 
describing a limit cycle oscillation, K is the coupling strength, h{ui) represents the 
feedback, and e is a unit vector with only one nonzero component: we have assumed 
that the feedback is additively applied to the system. 

When the coupling is weak, by treating the second term as a small perturbation 
the system is reduced to the phase model ([1]) [21]. In this phase description ([1]), 
synchronization behavior depends solely on the coupling function r(0), and therefore 
one can control the synchronization behavior of the system if the coupling function is 
freely given. It has been shown [191 EH] that this can be done by applying a properly 
designed external feedback signal h{ui) to the oscillators system. This method relies 
on the fact that the coupling function is the convolution of the feedback h and the 
phase response function Z{(f)), which characterizes the sensitivity of the phase to a weak 
external perturbation (see Appendix A ). 

Since oscillatory media can be regarded as a population of oscillators that are 
diffusively connected, we argue that the same method works for shaping the coupling 
function in the phase description of oscillatory media. Consider a (i-dimensional 
react ion- diffusion system with a global coupling: 

K f 

dtu = F{u) + DV'^u + / h{u)dx, (3) 

where u{x, t) is the state vector, D denotes the diffusion matrix, F{u) is a reaction term 
that generates a limit cycle oscillation with the frequency cu, K is the coupling strength, 
e is the same as above, and h{u) represents the feedback integrated over the entire 
space S. We assume that the system is Benjamin- Feir stable, i.e., the system undergoes 
spatially uniform oscillation when K = Q (external control is absent). Following the 
standard procedure (see Appendix A), we obtain 

dt(t){x, t) = uj + aV^cj) + (3{V(t)f + jj T{(P{x) - (t){x'))dx', (4) 

where 4>{x, t) is the phase of local oscillation; a, and (3 are constants determined by the 
property of the oscillatory medium. As in the case of discrete oscillators, the coupling 
function r(0) can be arbitrarily shaped by properly designing the feedback signal h. 

For discrete oscillators, ri-cluster states can be generated from the coupling function 
that contains v}^ harmonics [22]. Even in the case of oscillatory media, the coupling 
function is expected to work in the same way to stabilize the clustering pattern. A 
distinct problem here, however, is that the spatial patterns are not solely determined 
from the global coupling but from the interplay between the diffusive coupling and 
the global coupling, which makes the analysis much more complicated. For example. 
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when clustering pattern forms, interfaces appear between the clusters due to the diffusive 
coupling. Then controlling of the interface motion is required to obtain desired clustering 
patterns. 

Hence we take the following strategy. We start from a coupling function with which 
discrete oscillators described by equation ([T]) exhibit clustering or desynchronization. We 
then study the phase diffusion equation (4) with this coupling function to find resulting 
spatial patterns. Once the relation between the pattern and the coupling function 
is obtained, the corresponding feedback function can be found by following the same 
procedure as the discrete oscillators case. The key to carrying out this strategy is to 
find proper couping functions that allows for analytical treatment of the phase equation. 
Although analytical treatment is easy for a simple coupling function as r({/)) = sin (0+6') 
with a parameter 6 [1], only poor spatial patterns appear with such a coupling function; 
higher harmonics in the coupling function are responsible for the formation of complex 
spatial patterns, including phase clustering behavior. In the next section we propose 
such analytically tractable couphng functions that produces clustering states and the 
desynchronized state. 



3. Analysis of the phase model 

In this section, we study a one-dimensional phase diffusion equation with a global 
coupling (jll). Here we propose coupling functions that yield interesting spatial patterns; 
We are especially interested in two cluster states with specified population ratios, equally 
populated multi-cluster states, and the desynchronized state. 



3.1. 2-cluster states: numerical investigation 

Here we focus on the 2-cluster states with an arbitrary population ratio. In particular, 
we look for well-defined 2-clusters, with the phases maximally separated by vr. 

Discrete oscillators are known to form various clustering states, the behavior entirely 
governed by the form of the coupling function [22]. In special, the following coupling 
function yields 2-cluster states with the two phase difference equal to tt: 

r(0) = sin0-7{sin(20 + ^) -sin^}, (5) 



where 7 and 6 are parameters (see Appendix B). Equation ([T]) with this coupling 



function has a family of 2-cluster solution with different population ratios of the two 
clusters, which are stable for some range of population ratios. Thus, starting from a 
random initial condition the system converges to a 2-cluster state with its population 
ratio determined from the initial condition. 

Using this coupling function, we numerically investigate the phase diffusion 
equation (jlj) in one dimension with the system size S = L, taking 6* as a control 
parameter and 7 = 0.3. This equation is solved with the fiux-free boundary condition by 
using the second-order Euler scheme, with spatial and time interval being set to Ax = 0.1 
and At = 0.01, respectively. We set L = 100, K = 0.1, a = 0.384 x 10"^. We assign 
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Figure 1. Phase diagram for 2-cluster states and tlie desynchronized state in the 
phase model with ([5]) for each d = (3/ a. The points are numerical data, and the lines 
are given analytically. Recurrent 2-cluster has two ways of destabilization: the phase- 
advanced cluster becomes unstable (I), or the phase-retarded cluster becomes unstable 
(II). Inset figures are phase profiles with 6 = 2.83 and (a) 9 = 0.23, (b) 9 = 0.28, (c) 
9 = 1.58. In (b) the pattern is not stationary; here a snapshot of budding a new cluster 
is shown. 



several values to (3 to make a phase diagram below. Otherwise we set f3 = 1.089 x 10~^. 
This special choice of a and /5 is for later comparison to the Brusselator model. 

Figure [1] shows the phase diagram obtained by varying 6 for each 6 = j3/a with 
several values of (3 and fixed a. As expected, there exists a finite range of stationary 
2-cluster states. Note that, as opposed to discrete oscillators, here the population ratio 
between the two clusters is uniquely determined for fixed S and 6. Increase (decrease) in 6 
widens the phase-advanced (retarded) region. At some critical value of 6 the stationary 
state becomes unstable, leading to the recurrent 2-cluster state, where the following 
process occurs in a repeated way [see figure [2]^a)]: After a long transient of a quasi- 
stationary 2-cluster state, a new cluster sprouts out of the phase- advanced (retarded) 
cluster. Then the two interfaces propagate and one of the clusters disappear, the system 
returning to the 2-cluster state. Such dynamics have been reported in CO oxidation 
model [IE], although investigated only numerically. 

To characterize the patterns, we introduce the order parameters (/ = 1, 2, . . .): 

(^1 = 2 j dxe-'^'/'^'^l (6) 

For 2-cluster states, \ai \ indicates an approximate population disparity between the two 
clusters, and 1 — |(T2| the ratio of the interface width to the system size L. Note the two 
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different timescales in figure [21(b), each corresponding to the emergence of a new cluster 
and a slow drift of the interface. 

As 9 exceeds the threshold around ±7r/2, the recurrent 2-clusters turn into the 
desynchronized state [figure [T](c)], where a2 almost vanishes. 



3.2. 2- cluster states: analytical investigation 

Here we analytically investigate the 2-cluster states numerically found above. The 
analysis can be done by taking L oo limit. We derive analytical forms of ai and 
CT2 as functions of 9 and give the stability boundaries of the stationary 2-cluster shown 
in figure [H We move to a co-rotating frame so that the phase of the phase-retarded 
cluster is fixed to = 0. 

As we can see from the numerical result, the profile of a 2-cluster state can be 
decomposed into three regions: the phase- retarded cluster denoted by A (</> = 0), the 
phase-advanced cluster denoted by 5, and the interface. Also, from the numerical 
observation it is implied that the instability leading to the recurrent 2-clusters appears 
from the clustered region, while the interface remains stable. Hence in the analysis 
below we assume that the interface does not contribute to the stability. This separation 
of the regions becomes well-defined for large L. When the interface width is negligible 
compared to the system size, the two order parameters become real. In particular, in 
the steady state, we have (T2 = 1, so that Ui is the only relevant order parameter. 

Consider the dynamics of the cluster A. Contribution of the interface comes from 
the global coupling represented as the integral in equation (jl]). Since the interface 
width is vanishingly small, the interface region itself does not affect the dynamics. The 
remaining effect of the interface comes indirectly through the interface motion that varies 
the population ratio of the two clusters. However, since the timescale of the interface 
motion is 0(1/L), as shown below, the population ratio can be treated as constant. 
Thus in this limit the dynamics of the clusters is independent of the interface motion. 
When the interface can be negligible, equation (jlj) has a solution = for x & A 
and 4>{x) = TT for x G 5, where the population ratio is arbitrarily given. 

The stability analysis can be performed in the same way as the discrete oscillators 
(see Appendix B ). The only difference is the contribution from the diffusive coupling, 
which turns out to be negligible in the large L limit. Two modes of fluctuation occurs 
in the 2-cluster state: inter-cluster and intra-cluster fluctuation. Inter-cluster mode is 
a fluctuation of the phase between the clusters, with each cluster oscillating uniformly. 
The eigenvalue associated with this mode is given by Ainter = — 1 — 27 cos 6*. Thus by 
choosing I7I < ^ we can keep this mode stable. On the other hand, intra-cluster mode, a 
fluctuation within a cluster can be unstable. The eigenvalues associated with the cluster 
A and B with the wavenumber k are given by A^j^j.^^ = —ak"^ + 2p — 1 — 27 cos 9 and 
■^intra = — + 1 — 2p — 27 COS 9, respectively, where p, the population ratio, is the area 
fraction of the cluster A and is related to ai through 2p — 1 = ai. The negative sign of 
/c^-terms implies that the diffusion always works as stabilizing the inter-cluster modes; 
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Figure 2. (a) Space-time plot of a recurrent 2-cluster for a = 0.384 x 10~^, 
/3 = 1.089 X 10~^ and 6 = 0.28. Gray levels represent the spatial derivative of (f>, 
with the black lines indicating the location of the interfaces, (b) Corresponding time 
series of the amplitude of the order parameters ai and a2- (c)-(f) Snapshots of the 
phase profile, each corresponding to the arrows in (b). 



the most unstable mode is the one with the smallest (but finite) wavenumber. Taking 
the large L limit, this smallest wavenumber is vanishingly small, so that the /c^-terms 
can be dropped from the expression of the eigenvalue. 

Thus the diffusion does not affect the stability, while the stability depends on 
the population ratio. To obtain the analytical expression of the population ratio, let 
us consider the interface dynamics. The two clusters A and B are treated as the 
fixed boundaries of the interface. Since the inter-cluster mode is stable, the boundary 
conditions of the interface profile are given by </>(— oo) = and (f){oo) = tt, and a2 is 
replaced by the steady-state value, (72 = 1. Then equation (jlj) becomes 

dt(j) = dlcj) + 5{d^(t)f + cTsin0 - 7(sin(20 + 6)- sin^), (7) 

where we have defined a = ai and 5 = and rescaled time and space as Kt — >■ t 
and y/Kj 

Suppose that there is a traveling solution = f{x — ct), where the interface velocity 
c is written as c = f owing to the fact that interface motion varies the population 
ratio. Multiplying ([7]) by dr^f and integrating over the entire space yields 
da 4 f 5 /",_ , 71^ sin 9) 

Therefore a has a stable solution formally written as 



a 



^ f , 7r7sin6' . , 

-^Jidjfdx-^^. (9) 

From ([8]) it is seen that the interface slowly moves with the timescale of 0(1 /L) toward 
this stable state. Thus the interface dynamics, or the time evolution of a, is decoupled 
from the rest. 

An explicit expression of the steady state solution of a can be obtained through 
perturbation expansion. First we seek for a stationary solution of ([7]). When 6 satisfies 
tan^ = —6, there exists an exact solution connecting (f){—oo) = and 0(oo) = vr: 

</)o(x) = 2arctane'''', (10) 
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Figure 3. (a)Relation between the real part of a and for a = 0.384 x 10~^ and 
/3 = 1.089 X 10~^. The open circles are numerical data; the dotted lines represent 
the boundary of the existence of stationary 2-cluster states, given by (fT^ : the solid 
line is given by ((HI); crosses a.t 9 = —1.40 and 9 = 0.23 are numerical data from 
the Brusselator model with corresponding parameters, (b) Relation between critical 
values of a and 9, denoted by a* and 9*. The filled and open circles correspond to the 
instabilities of phase-advanced and phase-retarded clusters, respectively, with different 
values of f3. The dotted line is p^ . 



where k = a/27 cos 6. It is easily verified from (Q tliat tliis solution satisfies a = 0. 
Then perturbation expansion can be performed in terms of a up to the first order. We 
obtain (see [Appendix CD : 

27 cos 9 



a 



{S + tanO), (11) 



where x(^) = ^"i+t/a coth7r(5. 

Substituting the expression of a into the intra-cluster eigenvalues, we obtain the 
stability condition of 2-cluster states. In the large L limit, the /c^-term in the eigenvalues 
vanishes and the stability boundary is given by 

a ± 27 cos ^ = 0. (12) 

Figure [3] shows the dependence of a on 6. We have plotted only the real part of a, while 
in our simulation the imaginary part is 0{K) and is negligible. Both equations ffTTl) and 
( IT2l) fit well with the numerical data. Moreover, substituting (fTTI) into (fT2|) yields 

tanO = -5tx{S), (13) 

which gives the threshold value of 6 for the stability of stationary 2-cluster states. Note 
that the stability condition is independent of 7. The theoretical lines given by ( fT3l) are 
in excellent agreement with numerical data in figure [H 

Thus, within the range of parameter 6 determined from (|T3l) . we can control a as 
a function of 6 via (fTTj) . 

Now the interpretation of the recurrent 2-cluster states are given as follows: given 
an initial condition, the system converges to a 2-cluster state with a slowly-moving 
interface. When a, varies through the interface motion, exceeds the threshold given by 
equation ( fT2|) . the intra-cluster mode becomes unstable and one of the clusters collapses. 
Since the inter-cluster mode remains stable, the system returns to a 2-cluster state with 
reduced a, the whole process repeated ad infinitum. 
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3.3. Desynchronized state 

We give a theoretical analysis of equation ([7]) for the desynchronized state. The perfect 
desynchronized state is defined such that all the order parameters vanish. However, in 
practice some order parameters remain finite because of the boundary effect for flux-free 
boundary conditions. 

Firstly, as an ideal case, let us assume the periodic boundary condition. Then the 
perfect desynchronized state can be given by = 2'kx/L. Linear stability analysis 
for this profile shows that for each mode with the wavenumber ki = 2ttI/L (/ > 1) the 
corresponding eigenvalue is A^g^yj^^, = Aq^ — akf — 2(3ikiki, where Aq^'' = — Aq^^ = 76*^, 
and Aq"^"* = 0. Hence the I = 2 mode loses the stability at 6' = ±| in the large L limit, 
which fits well with numerical data in figure [H Note that, if the diffusive coupling is 
absent, only / = 1 and / = 2 modes are stable, / > 3 modes being neutral. 

Since the boundary is not periodic but fiux-free in the present case, the profile 
deviates from the linear one, as shown in figure [H^c), Accordingly, the steady state 
values of ai are shifted from zero by 0(1/ L) (order of the width of the boundaries). 
Note that in the linear regime the main contribution to ai comes from the mode ki. The 
modes / = 1 and / = 2 have the eigenvalues of order 0(1) as seen above, and thus ai and 
(72 remains to be 0(1/L). On the other hand, since / > 3 modes have the eigenvalues of 
0(1/L^), nonlinear effects of order 0(1/L^) coming from the terms such as (Tia2 makes 
/ > 3 modes grow up to 0(1). Therefore, in order to get better desynchronized state, 
we need to add as many higher harmonics as possible, as demonstrated in section HI 



3.4. Multi- cluster states 

The arguments of 2-cluster states can be extended to n-clusters in the following way. 
Consider the following coupling function: 

r(0) = sin m0 — 7{sin(n0 -\- 9) — sin 9}. (14) 

m=l 

This coupling function, when introduced to discrete oscillators, creates stable equally- 



populated n-clusters with the phases evenly separated ( [Appendix Bj ). Let us find a 
stationary, equally populated n-cluster solution of (jl]) with f|T4|) . Such a solution satisfies 
o'm = (m < n) and (T„ = 1, and hence only v}^ harmonic remains in (j4]). Then, by 
choosing 9 so as to satisfy 

5 = tan^, (15) 

we have a solution with each cluster separated by ^ and all the n — 1 interfaces having 
the same interface profile given by 0(a:) = - arctanexp(K„x) with k„ = ^/n^ cos 9. This 



state is stable against inter- and mira-cluster fluctuations (see Appendix B). 

For the stability of the desynchronized state, the same argument as in the above 
n = 2 case holds and the stability boundary is given by 6* = ±7r/2. 
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4. Numerical confirmation with the Brusselator model 

The above analytical expressions are used for the control of oscillatory media. As a 
model system of oscillatory media, we adopt the Brusselator model: 
du K f 

— = DuV\ + A-{B + l)u + u\ + j^ h{u, v)dx, (16) 

^ = D„V^v + Bu- v?v. (17) 
at 

The parameters A, i?, and are chosen in such a way that the system exhibits 
stable uniform oscillation; we set A = 1.6, B = 5.0, = 0.01, and = 0. The 
corresponding parameters in the phase model are a = 0.384 x 10~^ and jS = 1.089 x 10~^. 
For the precision that assures the validity of the phase description, we set K = 0.001 
and L = 1000 (equivalent to K = 0.1 and L = 100 in the phase model). Note that 
while for convenience of numerical simulation we have set i^t, = 0, we may also consider 
nonzero Dy, which simply results in the variation in the values of a and /3. 
As a feedback function we propose the following: 

M 

h{u, v) = h{(j)) = kn cos{n(j){u, v) — ipn) 5 (18) 

71=0 

where 0(m, v) is the phase of the limit cycle oscillation 0, and the parameters kn and 
ipn are the feedback intensity and the phase shift of the n^^ feedback term, respectively. 
The coupling function is obtained from h{(j)) and the phase response function Z{(j)), 
which characterizes the sensitivity of the phase to a weak external perturbation (see 



Appendix A ). By expanding Z{(j)) = zi cos{n(j) + xi) ^ the coupling function is written 



as 

r(0) = EV'°'^^'^ + ^' + ^'^- ^^^^ 

1=0 

In principle, as long as zi is finite, we can assign any value to l^^ harmonics of the 
coupling function by choosing appropriate values for ki and ipi. The advantage of 
using ( fTSj) is that the relation between parameters in the coupling function and the 
feedback parameters kn,ipn is given in a simple manner. (We could also use as the 
feedback h{u, v) a polynomial of u with multiple time delays [IHIEO], but in that case the 
relation is represented as a nonlinear function and the parameters need to be calculated 
numerically.) Hence, given a coupling function, we can calculate the corresponding 
feedback parameters by measuring phase response function Z{(j)). Table [1] shows Z{(j)) 
and the feedback parameters corresponding to (fT4l) for n = 2 and n = 5. In the following 
numerical investigation we set 7 = 0.3. 

I In our simulation, the phase 4){u, v) can be obtained directly from u and v in the following way: We 
first define the phase on the unperturbed {K = 0) limit cycle so that the phase evolves with a constant 
velocity, which can be done numerically. We then define the phase of a point(u, v) off the limit cycle 
by the phase of the nearest point on the limit cycle. 
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Table 1. Numerically obtained response function Z{(f)) ~ J^i cos{l(t) + xi) for the 
Brusselator with A ~ 1.6 and B = 5.0, and the feedback parameters {fc;} and {ijji} in 
(|18p producing the coupling functions given by (|14p with n ~ 2 and rt = 5. For I > n, 
ki and "i/); are equal to zero. 
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Figure 4. Two dimensional stationary pattern (uniform rotation subtracted) of the 
Brusselator model with global feedback containing 5 harmonics: (a) nearly equally 
populated 5-cluster {6 = —0.848) and (b) desynchronized state {6 = —1.580). (c) 
Section plots of (a) (top) and (b) (bottom), each from bottom left to top right. 

First we study 2-cluster states in the one- dimensional case with the feedback 
parameters corresponding to n = 2. We have confirmed that, for several parameter 
values of 6 we can observe stationary 2-clusters, recurrent 2-clusters, and desynchronized 
states, with the order parameter values predicted by the phase model (deviation of order 
O(10~^)). As an example, in figure [3](a), numerically obtained critical values of the real 
part of (Ti (denoted by "+") are superimposed on the data from the phase model, which 
are in good agreement with the corresponding phase model, with deviations O(10~^). 

Next, we use ( fT^ forn = 5 to produce the equally populated 5-cluster state and the 
desynchronized state in the two-dimensional case. The 5-cluster is shown in figure |l](a), 
the parameter 6 given by ( |T5i) with n = 5. The order parameters are {a^l = 0.739 
and Icr^l ~ O(10~^) for / < 5, indicating that the five clusters are well-defined and 
approximately equally populated. In figure 111(b) we have a desynchronized state with 6 
just below the threshold {6 = —1.58), where \ai\ ~ O(10~^) for / < 5 and O(10~^) for 
/ > 5. Note that the degree of desynchronization becomes better than the one for n = 2 
shown in section [3l we can make a better desynchronized state by adding appropriate 
higher harmonics. 
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5. Discussion: experimental realization of the theory 

To apply our method to experimental systems, we need to find the constants a and /3, 
and the response function Z{(f)). Since it is generally expected that the target pattern 
appears in the oscillatory media (due to inhomogeneities, or by applying a manual 
stimulus) [21], P can be measured by using the target pattern, assuming its phase 
profile as (f){x,t) = Qt + k\x\ (the origin is on the center of the target pattern) with 
the measurable quantities k and Q: Substituting this expression into (jlj) with K = 
we obtain (3 = [Q — uj)/k'^. The decay rate of the local perturbation from a uniform 
oscillation gives a. The response function 2'(0) can be measured by perturbing the 
system through a global parameter, to which we also apply the global feedback. To use 
( IT8|) . instantaneous measurement of the phase (j){x,t) at each spatial point is needed. 
If at least one quantity of an oscillator is observable, this can be done by, for example, 
constructing a delayed coordinate. 

Moreover, the following things should be taken into account for experimental 
realization of our theory. First, feedback must be weak for the precision of the phase 
description. This implies that the system size should be large enough to obtain well- 
defined cluster states even under weak feedback: the width of the interface is 0{\J D/ K), 
which must be sufficiently smaller than the linear dimension L. Also, to make a coupling 
function containing large enough higher harmonics with weak feedback, oscillation is 
better to be relaxation type: then the response function has higher harmonics with 
large amplitudes, so that we can keep the feedback signal weak to realize a desired 
coupling function [see the expression of r(</)) f|T9|) ]. Second, the emergence of phase 
singularity leads to the breakdown of the phase description and must be avoided. 

We have checked in our preliminary numerical simulations that the multi-cluster 
states and the desynchronized state are robust against noise. Thus we are convinced 
that our proposed method works in experimental systems. 

6. Concluding remarks 

We have proposed a theoretical framework for designing spatial patterns in oscillatory 
media. When a certain pattern is found in a phase model with a specific coupling 
function, the same pattern can be realized in oscillatory media by applying a properly 
constructed nonlinear feedback. In this paper, we found analytically tractable coupling 
functions that enables us to quantitatively control the spatial patterns. Using these 
coupling functions, we investigated the phase equation with the global coupling and 
found the parameter regions where the following patterns stably exist: 2-cluster 
states with specified population ratios, equally populated multi-cluster states, and the 
desynchronized state. In the case of 2-clusters, we gave analytical expression of the 
population ratio of the two clusters as the function of a feedback parameter. We also 
proposed a simple form of the nonlinear feedback function to make the calculation of 
the feedback parameters easier. We exemplified all these results using the Brusselator 
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model and succeeded to reproduce the patterns predicted by the phase model. Since 
our method is based on the measurable quantities only, it is expected that the method 
is verified in a real experiment. 

The desynchronized state deserves further remark. Our results show that even 
in oscillatory media one can drive the system into the desynchronized state, as well 
as in discrete oscillators [19]. Such a control is not only of medical [23] , but also 
potentially of industrial interest; for example, it would be beneficial when constant 
output from oscillatory catalytic reaction is desirable. 

Further investigation of the phase model with other coupling functions is of great 
interest for controlling more complex patterns, although our method is limited to 
oscillatory system and cannot be applied to some typical spatial patterns such as 
the Turing pattern. Also, investigating the control of Benjamin-Feir unstable systems 
by replacing the phase diffusion equation with Kuramoto-Sivashinsky equation will be 
interesting both in a theoretical sense and for application. 
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Appendix A. Derivation of the phase model 

In this Appendix we derive the phase model (jlj) from a reaction-diffusion system with a 
global feedback. The system is assumed to undergo spatially uniform oscillation when 
external control is absent (namely, the system is Benjamin-Feir stable [21]). Dynamical 
evolution of a c?- dimensional oscillatory medium is described by a reaction-diffusion 
equation: 



Note that here we consider a general situation, where the global feedback is introduced 
through a global parameter q. In equation (I3l), and in References [19], [20], the feedback 
is simply applied additively. External feedback is applied to q as 



where go and K > Q are constants. By assumption, dtu = F{u] go) yields a limit-cycle 
oscillation, with its solution denoted by n = Wo(^)- The function p{t) describes a global 
feedback signal, given by 



where h{u) is some feedback function. The integration is taken over the entire space 
and S is the volume of the system. (Various functions can be considered for h. Our 
particular choice has been given in equation [181 ) 



dtU = F{u; q) + tv'^u. 



(A.l) 



q{t) = go + Kpit) 



(A.2) 




(A.3) 



Design principle of multi-cluster and desynchronized states in oscillatory media 



14 



As we have assumed, feedback intensity K is small, so that by dropping 0{K'^) 
equation flA.ip can be approximated by 

dtu = F{u; go) + DV'u + Kp{t)f{u), (A.4) 

where f{u) = {dF/dq)g=gg. When / is independent of u, the global parameter q 
appears additively and the system reduces to equation ([3]). 

When a spatial pattern emerges for small K > 0, the spatial variation, and thus 
V^it, is expected to be small, vanishing as K ^ 0. Thus, in addition to the feedback 
term, we may treat the diffusion term as small perturbations to the limit cycle (this is 
the case in our simulation, where the interface width is 0{\J D / K)), and therefore the 
diffusion term is the same order as the feedback). Then, following a standard method 
developed by Kuramoto [2T], we can derive a closed description for the phase variable 
for our oscillatory medium. As is usually adopted, the phase (t){u) is defined so as to 
satisfy (9^0 ■ F{u] go) = ^- Substituting this relation into the identity dtcj) = du4> ■ dtU, 
we obtain 

dt(l) = uj + d^cj) ■ jw^u + Kp{t)f{u)^ . (A.5) 

At the lowest order of we can replace u with the value on the limit cycle uq. Then 
the equation above is expressed only in terms of 0. After averaging ( 1A.5I) over one 
period of oscillation, we arrive at equation (jl]), where a, /3, and r(0) are written as 

^^l.jyzim% (A.6) 

^^^/"d^ZWD^, (A.7) 

r(0-0') = 7^/ dAZ(0 + A)/i(0' + A). (A.8) 

^TT Jo 

Here, the phase response function Z{(j)) = Z{(j)) ■ f{(f)), defined as the response to the 
global parameter g, and the "bare" response function Z{(j)) = du(p\u=uoy are evaluated 
on the unperturbed limit-cycle orbit. 

Expanding Z(0) and h{(j)) as Z(0) = TZ-oo ^i^^'*' and /i(0) = T.Z-M^ie'^'^ 
respectively, we obtain 

M 

r(0) = zih-ie''t (A.9) 

l=-M 

Hence the coupling function containing up to the M^^ harmonics can be generated by 
determining /i(0) up to the M**^ harmonics, as long as zi has a finite value p!9l [20] . 

Appendix B. 2-cluster states for the coupled oscillators 



Here we show that the collection of discrete oscillators interacting through the coupling 
function given by ([5]) can exhibit 2-cluster states with their phases separated by vr. 
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Appendix B.l. Steady-state 2- cluster solution 

Consider a set of TV identical oscillators with the frequency u. The dynamics is written 

as 

K ^ 
i=i 

To produce n-cluster states, it is sufficient that the coupling function contains up to the 
n^^ harmonics [22]; linear stability analysis shows that the harmonics smaller than n 
does not contribute to the stability of n-cluster states, and the n*^ harmonics works in a 
similar way to n:l periodic forcing [1]. In special, to observe 2-cluster states, one needs 
to prepare the coupling function such that the ffist harmonics destabilizes the 1-cluster, 
i.e., perfect synchronization, and the second assures the 2-cluster. Thus the coupling 
function for 2-cluster states can be written as 

r(0) = sin(0 + ^i) - 7 sin(2(/. + 62). (B.2) 

The amplitude of the ffist harmonics can be absorbed into the coupling constant K. 
Also, for later convenience we choose the negative sign for the second harmonics. 

Assume that the oscillators form a 2-cluster state, where A^^ oscillators belong 
to the cluster A with the phase (pA and A^ — Na to the cluster B with (pB- In the 
phase-locking state (0j = Vt for all z), we get 

fi = pr(o) + (i-p)r(^), (B.3) 

n = pT{-ij)+pT{{)), (B.4) 

where ip = (j)A — 4>b, P = Na/N and Vt = uj + sin^i — 7 sin ^2 is the frequency of the 
clusters. Then ip satisfies 

{2p - l)r(O) + (1 - p)V{^p) - pV{-^p) = 0. (B.5) 

When we choose 61 = 0, ( ]B.5|) has a solution ip = n ioi any p. In special, when I7I < |, 
ip = IT is the only solution except for ^^ = 0, the single cluster solution. If 9i 7^ 0, the 
phase difference is shifted from n except for p = |. 

Appendix B.2. Linear stability analysis 

We perform the linear stability analysis by expanding (pj as (pj = (p^^^ + where p^j'^ = 
for j G A and p^j'^ = ir for j G B. 

First we consider the mter-cluster mode, where the fluctuation is uniform in each 
cluster. In this case we can write ^j^A = and C,j£B = Cb- The mode C,a obeys 

U = K{l-p)T'in)i^A-^B), (B.6) 
^^ = KpT'i-7r)i^B-U)- (B.7) 

Changing the variables as = .^a ± C,b, and using r'(7r) = r'(— tt) = — 1 — 27 cos 6', we 

obtain 

e+ = -2p)(l + 27cos^)e-, (B.8) 
i_ = -K{1 + 27 cos 9)^^. (B.9) 
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The zero mode ^+ represents uniform rotation along with the hmit cycle. On the other 
hand, corresponds to the inter-cluster mode, which is stable regardless of 9 for I7I < i. 

Next we consider the mtra-cluster mode, where the fluctuation occurs within each 
cluster and the spatial average of ^ within each cluster is zero. We get 

OeA = K{2p - 1 - 27cos^)eieA, (B.IO) 
ij^B = Kil -2p- 2-fcos9)^jeB- (B.ll) 

(B.12) 

Thus the intra-cluster fluctuation the eigenvalues A^i = K{2p — 1 — 27 cos 6') and 
= A'(l — 2j9— 27 cos 9). These modes can be destabilized depending on the population 
ratio p. In special, when \9\ > 11/2, either or is positive for any p. 

Similarly, the coupling function ( iMl) produces n-cluster solutions 0j = 2nl/n 
(/ = 0, ■ ■ ■ , n — 1) for arbitrary population ratio, as can be checked by direct substitution. 
Stability is studied analogously with the n = 2 case above, the intra- and mter-cluster 
eigenvalues given by 



>^± = -K['- + n^cos9], (B.13) 



n 



2 

AS.a = (i - + ^^^°^^) ' 

where pm is the fraction of the m}^ cluster. When the clusters are equally populated, 
'^intra ~ —Kwy cos9 and is stable for \9\ < tx 12. Conversely, when 16*1 > 7r/2 at least one 
of n clusters has positive A^^^j.^^ and the ra-cluster state is no longer stable. 

Appendix C. Derivation of (1111) 

We expand the parameter 9 and the profile (^{x) in terms of a as follows: 

tan^ = -5 + (T/ii + 0((t2), (C.l) 
<\>{x) = M^) + ^M^) + 0{a^). (C.2) 

Substituting these expressions into ([7]) yields the following linearized equation: 

C(f)i{x) = —27 cos 6^ sin 00 (a^) — /Wi sin^ 0o(a^), (C.3) 

where the linearized operator C is given by 

C = K~^dl - cos2(j)o{x) + 26sm(f)o{x){K^^d^ - cosM^))- (C.4) 
The adjoint operator is written as 

£^ = — cos2</)o(x) — 25 sm(f)o{x){K~^dx + 2 cos</)o(x)). (C.5) 

It is verified by direct calculation that C'' has the zero eigenfunction 
exp[250o(a;)]sech/tx. The solvability condition reads 

dx exp [2(500 (x)]sechKx (— 27cos6'sin0o(a;) — /^i sin^ 0o(a;)) = 0, (C.6) 
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which yields 

_ J(l + ^^)coth7r^ 

7cose(l + 452) • ^^-'^ 

Using the relation tan6' = —5 + cr/ii, we obtain flTT]) . 
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